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ABSTRACT 

Dissipationless collapses in Modified Newtonian Dynamics (MOND) are studied by using a new 
particle- mesh N-body code based on our numerical MOND potential solver. We found that low 
surface-density end-products have shallower inner density profile, flatter radial velocity-dispersion 
profile, and more radially anisotropic orbital distribution than high surface-density end-products. 
The projected density profiles of the final virialized systems are well described by Sersic profiles with 
index m < 4, down to m ~ 2 for a deep-MOND collapse. Consistently with observations of elliptical 
galaxies, the MOND end-products, if interpreted in the context of Newtonian gravity, would appear 
to have little or no dark matter within the effective radius. However, we found impossible (under 
the assumption of constant mass-to-light ratio) to simultaneously place the resulting systems on the 
observed Kormendy, Faber- Jackson and Fundamental Plane relations of elliptical galaxies. Finally, the 
simulations provide strong evidence that phase mixing is less effective in MOND than in Newtonian 
gravity. 

Subject headings: gravitation — stellar dynamics — galaxies: kinematics and dynamics — galaxies: 
elliptical and lenticular, cD - galaxies: formation — methods: numerical 



1. INTRODUCTION 

In Bekenstein & Milgrom's (1984, hereafter BM) La- 
grangian formulation of Milgrom's (1983) Modified New- 
tonian Dynamics (MOND), the Poisson equation 



AnGp 



(1) 



for the Newtonian gravitational potential </>n is replaced 
by the field equation 



V ■ 



= AirGp, 



(2) 



where a,Q ~ 1.2 x 10 _10 ms -2 is a characteristic accel- 
eration, ||...|| is the standard Euclidean norm, <p is the 
MOND gravitational potential produced by the density 
distribution p, and in finite mass systems — * for 
|jx|j — ► oo. The MOND gravitational field g experienced 
by a test particle is 



g = -w, 

and the function p is such that 

ryj > [1 for y > 1; 
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(3) 



(4) 



throughout this paper we use 

Kv) = 



V 1 + y 2 

In the so-called 'deep MOND regime' (hereafter 
dMOND), describing low-acceleration systems (||V</>|| <C 
do), p(y) — y and so equation ^ simplifies to 



V • (||V0||V0) = 4irGa p. 



(6) 



The source term in equation ([2]) can be eliminated by 
using equation (fT]), giving 

|V0| 



X7<f> = V^n + S, 



(7) 



where S = curlh is a solenoidal field dependent on p 
and in general different from zero. When S = equa- 
tion {7} reduces to Milgrom's (1983) formulation and can 
be solved explicitly. Such reduction is possible for config- 
urations with spherical, cylindrical or planar symmetry, 
which are special cases of a more general family of strat- 
ifications (BM; Brada & Milgrom 1995). Though the 
solenoidal field S has been shown to be small for some 
configurations (Brada & Milgrom 1995; Ciotti, Londrillo 
& Nipoti 2006, hereafter CLN), neglecting it when simu- 
lating time-dependent dynamical processes has dramatic 
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effects such as non-conservation of total linear momen- 
tum (e.g. Felten 1984; see also Section [XT]) . 

Nowadays several astronomical observational data ap- 
pear consistent with the MOND hypothesis (see, e.g., 
Milgrom 2002; Sanders & McGaugh 2002). In addition, 
Bekenstein (2004) recently proposed a relativistic version 
of MOND (Tensor- Vector-Scalar theory, TeVeS), mak- 
ing it an interesting alternative to the cold dark matter 
paradigm. However, dynamical processes in MOND have 
been investigated very little so far, mainly due to difficul- 
ties posed by the non-linearity of equation ([2]) . Here we 
recall the spherically symmetric simulations (in which 

5 = 0) of gaseous collapse in MOND by Stachniewicz 

6 Kutschera (2005) and Nusser & Pointecouteau (2006). 
The only genuine three-dimensional MOND N-body sim- 
ulations (in which equation [2] is solved exactly) are those 
by Brada & Milgrom (1999, 2000), who studied the sta- 
bility of disk galaxies and the external field effect, and 
those of Tiret & Combes (2007). Other attempts to 
study MOND dynamical processes have been conducted 
using three-dimensional N-body codes by arbitrarily set- 
ting S = 0: Christodoulou (1991) investigated disk sta- 
bility, while Nusser (2002) and Knebe & Gibson (2004) 
explored cosmological structure formation 1 . 

In this paper we present results of N-body simula- 
tions of dissipationless collapse in MOND. The simula- 
tions were performed with an original three-dimensional 
particle-mesh N-body code, based on the numerical 
MOND potential solver presented in CLN, which solves 
equation ^ exactly. These numerical experiments arc 
interesting both from a purely dynamical point of view, 
allowing for the first time to explore the relaxation pro- 
cesses in MOND, and in the context of elliptical galaxy 
formation. In fact, the ability of dissipationless collapse 
at producing systems strikingly similar to real ellipti- 
cals is a remarkable success of Newtonian dynamics (e.g., 
van Albada 1982; Aguilar & Merritt 1990; Londrillo, 
Messina & Stiavelli 1991; Udry 1993; Trenti, Bertin & 
van Albada 2005; Nipoti, Londrillo, & Ciotti 2006, here- 
after NLC06), while there have been no indications so 
far that MOND can work as well in this respect. Here 
we study the structural and kinematical properties of 
the end-products of MOND simulations, and we com- 
pare them with the observed scaling relations of ellipti- 
cal galaxies: the Faber-Jackson (FJ) relation (Faber & 
Jackson 1976), the Kormendy (1977) relation, and the 
Fundamental Plane (FP) relation (Djorgovski & Davis 
1987, Dressier et al. 1987). 

The paper is organized as follows. The main features 
of the new N-body code are presented in Section [2J while 
Section [3] describes the set-up and the analysis of the 
numerical simulations. The results are presented in Sec- 

1 Cosmological N-body simulations in the context of a relativistic 
MOND theory such as TeVeS have not been performed so far. 



tion[4]and discussed in Section [5l 

2. THE N-BODY CODE 

While most N-body codes for simulations in New- 
tonian gravity are based on the gridless multipole ex- 
pansion treecode scheme (Barnes & Hut 1986; see also 
Dehnen 2002), the non-linearity of the MOND field equa- 
tion © forces one to resort to other methods, such as the 
particle- mesh technique (see Hockney & Eastwood 1988). 
In this approach, particles are moved under the action of 
a gravitational field which is computed on a grid, with 
particle-mesh interpolation providing the link between 
the two representations. In our MOND particle-mesh 
N-body code, we adopt a spherical grid of coordinates 
(r, $, cp), made of N r x N# x N v points, on which the 
MOND field equation is solved as in CLN. Particle-mesh 
interpolations are obtained with a quadratic spline in 
each coordinate, while time stepping is given by a classi- 
cal leap-frog scheme (Hockney & Eastwood 1988). The 
time-step At is the same for all particles and is allowed 
to vary adaptively in time. In particular, according to 
the stability criterion for the leap-frog time integration, 
we adopt At = ^/ymax |V 2 0|, where r\ < 0.3 is a dimcn- 
sionless parameter. We found that 77 = 0.1 assures good 
conservation of the total energy in the Newtonian cases 
(see Section l3"7Tj) . In the present version of the code, all 
the computations on the particles and the particle-mesh 
interpolations can be split among different processors, 
while the computations relative to the potential solver 
are not performed in parallel. The solution of equa- 
tion @ over the grid is then the bottleneck of the sim- 
ulations: however, the iterative procedure on which the 
potential solver is based (see CLN) allows to adopt as 
seed solution at each time step the potential previously 
determined. 

The MOND potential solver can also solve the Poisson 
equation (obtained by imposing fi = 1 in equation^, so 
Newtonian simulations can be run with the same code. 
We exploited this property to test the code by running 
several Newtonian simulations of both equilibrium distri- 
butions and collapses, comparing the results with those 
of simulations (starting from the same initial conditions) 
performed with the FVFPS treecode (Londrillo, Nipoti 
& Ciotti 2003; Nipoti, Londrillo & Ciotti 2003). One of 
these tests is described in Section l4.1.2l 

We also verified that the code reproduces the New- 
tonian and MOND conservation laws (see Section l3~Tj) : 
note that the conservation laws in MOND present some 
peculiarities with respect to the Newtonian case, so we 
give here a brief discussion of the subject. As already 
stressed by BM, equation ^ is obtained from a vari- 
ational principle applied to a Lagrangian with all the 
required symmetries, so energy, linear and angular mo- 
mentum are conserved. Unfortunately, as also shown by 
BM, the total energy diverges even for finite mass sys- 
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TABLE 1 

Time, velocity, and energy units for Newtonian and 
MOND (subscript n), and dMOND (subscript d) N-body 
simulations. 

Un = rl /2 {GM,)' 1 / 2 t, d =r,(GM,tt )- 1 / 4 
v* n = (GM^/V; 1 ' 2 u» d = (GA-Uao) 1 / 4 



E tn = GMt, 



terns, thus posing a computational challenge to code val- 
idation. We solved this problem by checking the volume- 
limited energy balance equation 



dtJvo 



8irG 



4ttG 



dV 



8^ 

Li— < V0, n > da, 
at 



(8) 



which is derived in Appendix [X] In equation §E§ Vq is 
an arbitrary (but fixed) volume enclosing all the system 
mass, k is the kinetic energy per unit volume, and 

f(y)E2f^M, (9) 
Jyo 

where yo is an arbitrary constant; note that only finite 
quantities are involved. Another important relation be- 
tween global quantities for a system at equilibrium (in 
MOND as in Newtonian gravity) is the virial theorem 



2K + W = 0, 



(10) 



where K is the total kinetic energy and W — Tr Wij is 
the trace of the Chandrasekhar potential energy tensor 



W ir , 



p(x)xi 



dx 



(11) 



(e.g., Binney & Tremaine 1987). Note that in MOND 
K + W is not the total energy, and is not conserved. 
However, W is conserved in the limit of dMOND, be- 
ing W — —(2/3)y/GaoM% for all systems of finite total 
mass M» (see Appendix [5] for the proof). As a con- 
sequence, in dMOND the virial theorem writes simply 
<7y = 4GM*ao/9, where ay = y/ 2K/ M* is the system 
virial velocity dispersion (this relation was proved for 
dMOND spherical systems by Gerhard & Spergel 1992; 
see also Milgrom 1984). In our simulations we also tested 
that equation (JTUJ) is satisfied at equilibrium, and that W 
is conserved in the dMOND case (see Sections 13.11 and 

3. NUMERICAL SIMULATIONS 

The choice of appropriate scaling physical units is an 
important aspect of N-body simulations. This is espe- 
cially true in the present case, in which we want to com- 
pare MOND and Newtonian simulations having the same 
initial conditions. As well known, due to the scale- free 
nature of Newtonian gravity, a Newtonian ./V-body sim- 
ulation starting from a given initial condition describes 



in practice oo 2 systems of arbitrary mass and size. Each 
of them is obtained by assigning specific values to the 
length and mass units, r* and M*, in which the initial 
conditions are expressed. Also dMOND gravity is scale 
free, because ao appears only as a multiplicative factor 
in equation ([6]), and so a simulation in dMOND gravity 
represents systems with arbitrary mass and size (though 
in principle the results apply only to systems with accel- 
erations much smaller than a ). MOND simulations can 
also be rescaled, but, due to the presence of the charac- 
teristic acceleration ao in the non-linear function /it, each 
simulation describes only oo 1 systems, because r* and 
M* cannot be chosen independently of each other. 

On the basis of the above discussion, we fix the physical 
units as follows (see Appendix [C] for a detailed descrip- 
tion of the scaling procedure). Let the initial density 
distribution be characterized by a total mass M* and a 
characteristic radius r* . We rescale the field equations so 
that the dimensionless source term is the same in New- 
tonian, MOND and dMOND simulations. We also re- 
quire that the Second Law of Dynamics, when cast in 
dimensionless form, is independent of the specific force 
law considered, and this leads to fix the time unit. As a 
result, Newtonian and MOND simulations have the same 
time unit t* a = r^ 2 (GA/*)~ 1 / 2 , while the natural time 
unit in dMOND simulations is £*d = r„(GM t ,ao)^ 1 ^ 4 . 
Note that MOND simulations are characterized by the 
dimensionless parameter k = GM*/r 2 ao, and scaling of 
a specific simulation is allowed provided the value of k is 
maintained constant. So, simulations with lower k val- 
ues describe lower surface-density, weaker acceleration 
systems; dMOND simulations represent the limit case 
k 1 , while Newtonian ones describe the regime with 
k» 1. With the time units fixed, the corresponding ve- 
locity and energy units are w* n = r*/t* n; v *d = r*/t*<n 
E m = M*uJ n , and E* d = Af*v 2 d (see Table 1 for a sum- 
mary). 

3.1. Initial conditions and analysis of the simulations 

We performed a set of five dissipationless-collapse N- 
body simulations, starting from the same phase-space 
configuration: the initial particle distribution follows the 
Plummer (1911) spherically symmetric density distribu- 
tion 

Pir) - ™: r l^ (12) 



47r(r 



2)5/2 



where is the total mass and r* a characteristic ra- 
dius. The choice of a Plummer sphere as initial con- 
dition is quite artificial, and not necessarily the most 
realistic to reproduce initial conditions in the cosmolog- 
ical context (e.g., Gunn & Gott 1972). We adopt such a 
distribution to adhere to other papers dealing with col- 
lisionless collapse (e.g., Londrillo et al. 1991; NLC06; 
see also Section [5l in which we present the results of a 
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Fig. 1.— Time evolution of 2K/\W\, K, W, and K + W for simulations D, Ml, and N. K, W, and K 
column), and _E* n (central and right columns). For clarity, the time axes are zoomed-in between and 10. 



W are in units of (left 



set of simulations starting from different initial condi- 
tions). The particles are at rest, so the initial virial ratio 
2K/\W\ — 0. What is different in each simulation is 
the adopted gravitational potential, which is Newtonian 
in simulation N, dMOND in simulation D, and MOND 
with acceleration ratio k in simulations Mk (k=1, 2, 4). 
For each simulation we define the dynamical time idyn 
as the time at which the virial ratio 2K/\W\ reaches its 
maximum value. In particular, we find td- 



vn 



2£*d in 

simulation D, and idyn ~ 2i* n in simulations N, Ml, M2 
and M4. We note that t dyn ~ GM* 5/2 (2\K + W\)~^ 2 in 
simulation N. 

Following NLC06, the particles are spatially dis- 
tributed according to equation (fT2|) and then randomly 
shifted in position (up to r*/5 in modulus). This arti- 
ficial, small-scale "noise" is introduced to enhance the 
phase mixing at the beginning of the collapse, because 
the numerical noise is small, and the velocity dispersion 
is zero (see also Section 4.2). As such, these fluctuations 
are not intended to reproduce any physical dumpiness. 

All the simulations (realized with N = 10 6 particles, 



and a grid with N r = 64, N$ — 16 and N v = 32) are 
evolved up t = 150id yn - In all cases the modulus of the 
center of mass position oscillates around zero with r.m.s 
< 0.1r»; similarly, the modulus of the total angular mo- 
mentum oscillates around zero 2 with r.m.s. <0.02, in 
units of r*M*v m (simulations Mk and N) and of ?'»M»w*d 
(simulation D). K + W in the Newtonian simulation and 
W in the dMOND simulation are conserved to within 2% 
and 0.6%, respectively. The volume-limited energy bal- 
ance equation §5§ is conserved with an accuracy of 1% in 
MOND simulations, independently of the adopted Vq. To 
estimate possible numerical effects, we reran one of the 
MOND collapse simulations (Ml) using N — 2 x 10 6 , 
N r = 80, N# = 24, and N v = 48: we found that the 
end-products of these two simulations do not differ sig- 



2 As an experiment we also ran a simulation, with the same 
initial conditions and parameter re as Ml, in which the force was 
calculated from equation J7) imposing S = 0. In this simulation 
the linear and angular momentum are strongly not conserved: for 
instance, the center of mass is already displaced by ~ 7r* after 

~ 30tdyn- 
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nificantly, as far as the properties relevant to the present 
work are concerned. 

The intrinsic and projected properties of the collapse 
end-products are determined as in NLC06. In partic- 
ular, the position of the center of the system is deter- 
mined using the iterative technique described by Power 
et al. (2003). Following Nipoti et al. (2002), we measure 
the axis ratios c/a and b/a of the inertia ellipsoid (where 
a, b and c are the major, intermediate and minor axis) of 
the final density distributions, their angle-averaged pro- 
file and half-mass radius rt- We fitted the final angle- 
averaged density profiles with the 7-model (Dehnen 1993; 
Tremaine et al. 1994) 



p(r) 



rT(r c + r) 4 -f ' 



(13) 



where the inner slope 7 and the break radius r c are 
free parameters, and the reference density po is fixed 
by the total mass M*. The fitting radial range is 
0.06 < r/r^ < 10. In order to estimate the importance of 
projection effects, for each end-product we consider three 
orthogonal projections along the principal axes of the in- 
ertia tensor, measuring the ellipticity e = 1 — b c /a c , the 
circularized projected density profile and the circularized 
effective radius R c = \/a c b c (where a e and b e are the ma- 
jor and minor semi-axis of the effective isodensity ellipse). 
We fit (over the radial range 0.1 < R/R c ^ 10) the cir- 
cularized projected density profiles of the end-products 
with the R}l m Sersic (1968) law: 



I(R) = I c cxp <^ -b{m) 



(14) 



where I c = I{R C ) and b(m) ~ 2m- l/3 + 4/405m (Ciotti 
& Bertin 1999). In the fitting procedure m is the only free 
parameter, because R c and I c are determined by their 
measured values obtained by particle count. In addition, 
we measure the central velocity dispersion <7oi obtained 
by averaging the projected velocity dispersion over the 
circularized surface density profile within an aperture of 
J? e /8. Some of these structural parameters are reported 
in Table 2 for the five simulations described above, as 
well as for three additional simulations, which start from 
different initial conditions (see Section [5]). 

4. RESULTS 

In Newtonian gravity, collisionless systems reach viri- 
alization through violent relaxation in few dynamical 
times, as predicted by the theory (Lynden-Bell 1967) 
and confirmed by numerical simulations (e.g. van Al- 
bada 1982). On the other hand, due to the non linearity 
of the theory and the lack of numerical simulations, the 
details of relaxation processes and virialization in MOND 
are much less known. Thus, before discussing the spe- 
cific properties of the collapse end-products we present a 



general overview of the time evolution of the virial quan- 
tities in our simulations, postponing to Section r4.2l a more 
detailed description of the phase-space evolution. In par- 
ticular, in Fig. [1] we show the time evolution of 2K/\W\, 
K, W, and K + W for simulations D, Ml, and N. In the 
diagrams time is normalized to idyn, so plots referring to 
different simulations are directly comparable (the values 
of tdyn in time units for the five simulations are given 
in Section 13. ip . In simulation N (right column) we find 
the well known behavior of Newtonian dissipationless col- 
lapses: 2K/\W\ has a peak, then oscillates, and eventu- 
ally converges to the equilibrium value 2K/\W\ = 1; the 
total energy K + W is nicely conserved during the col- 
lapse, though it presents a secular drift, a well known 
feature of time integration in N-body codes. The time 
evolution of the same quantities is significantly different 
in a dMOND simulation (left column). In particular, the 
virial ratio 2K/\W\ quickly becomes close to one, but is 
still oscillating at very late times because of the oscilla- 
tions of K, while W is constant as expected. As we show 
in Section [4~2l these oscillations are related to a peculiar 
behavior of the system in phase space. Finally, simula- 
tion Ml (central column) represents an intermediate case 
between models N and D: the system starts as dMOND, 
but soon its core becomes concentrated enough to enter 
the Newtonian regime. After the initial phases of the 
collapse, Newtonian gravity acts effectively in damping 
the oscillations of the virial ratio. Overall, it is apparent 
how the system is in a "mixed" state, neither Newto- 
nian (K + W is not conserved) nor dMOND (W is not 
constant). 

4.1. Properties of the collapse end-products 
4.1.1. Spatial and projected density profiles 

We found that all the simulated systems, once viri- 
alized, are not spherically symmetric. However, while 
the dMOND collapse end-product is triaxial (c/a ~ 0.2, 
b/a ~ 0.4), MOND and Newtonian end-products are 
oblate (c/a ~ c/b ~ 0.5). The ellipticity e of the pro- 
jected density distributions (measured for each of the 
principal projections) is found in the range 0.5 — 0.8 in 
D, and - 0.5 in Ml, M2, M4 and N. These values are 
consistent with those observed in real ellipticals, with the 
exception of et in model D (see Table 2), which would 
correspond - if taken at face value - to an E8 galaxy. 
These result could be just due to the procedure adopted 
to measure the ellipticity (see Section l3~Tj) , however we 
find interesting that dMOND gravity could be able to 
produce some system that would be unstable in New- 
tonian gravity. We remark that a similar result, in the 
different context of disk stability in MOND, has been 
obtained by Brada & Milgrom (1999). 

In order to describe the radial mass distribution of the 
final virialized systems, we fitted their angle-averaged 
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TABLE 2 
End-product properties. 





K 


c/a 


b/a 


7 


r c /r h 




m 6 


m c 


e a 


e& 


e c 


D 
Ml 

M2 
M4 


1 

2 
4 


0.21 
0.47 
0.48 
0.47 


0.41 
0.85 
0.92 
0.90 


u - 1 '-0.17 
1 24+ ' 40 

1 45+ 026 

1 54+0.26 


n 97 +u.uy 

u - z '-0.02 
44+°- 20 

0.531°;" 
0.581°;?° 


2.87 ±0.01 
3.20 ±0.07 
3.38 ±0.08 
3.55 ±0.10 


2.50 ± 0.03 
3.00 ±0.09 
3.24 ±0.08 
3.40 ±0.11 


2.16 ±0.02 
3.07 ±0.13 
3.28 ±0.12 
3.34 ±0.15 


0.48 
0.42 
0.49 
0.51 


0.80 
0.51 
0.45 
0.45 


0.58 
0.17 
0.08 
0.10 


N 




0.45 


0.91 


1 69+ - 13 
1 ' oa -0.15 


0.74lg;l| 


4.21 ±0.07 


4.35 ±0.08 


3.96 ±0.13 


0.48 


0.55 


0.12 



T7 - 025 045 0.721^ 0.37^^ 3.06 ± 0.06 2.90 ± 0.04 2.71 ±0.08 044 076 0~56 
M' 20 0.42 0.83 1.26t° ^ 0.47^^1 3.41 ± 0.09 3.36 ± 0.06 3.20 ±0.13 0.49 0.57 0.16 
N' - 0.45 0.93 1.78±";^ 0.781°;^ 4.29 ±0.10 4.56 ±0.15 4.19 ±0.22 0.51 0.55 0.09 



First column: name of the simulation, k = GM* /r 2 ao: acceleration ratio, c/a and b/a: minor-to-major and intermediate-to-major axis 
ratios. 7, r c : best-fit 7-model parameters. m a , rribi m c and e a , e^, e c : best-fit Sersic indices and cllipticities for projections along the 
principal axes. 




Fig. 2. — Angle-averaged density, radial velocity-dispersion and anisotropy-parameter profiles (from bottom to top) for the end-products 
of simulations D, Ml, and N. Dotted lines in the bottom panels represent p <x r _1 profiles, which are shown for reference. Empty circles 
in the right column show the corresponding profiles obtained with the FVFPS treecode from the same initial conditions. 
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Fig. 3. — Line-of-sight velocity-dispersion profiles (top), circularized projected density profiles and residuals of the Sersic fit (bottom) of 
the end-products of simulations D, Ml, and N (squares; 1-tr error bars are always smaller than the symbol size). The dotted lines are the 
best-fitting Sersic models. 
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density profiles with the 7-model (|13p over the radial 
range 0.06 < r/r n < 10. The best-fit 7 and r c for the 
final distribution of each simulation are reported in Ta- 
ble 2 together with their la uncertainties (calculated 
from A\ 2 = 2.30 contours in the space 7 — r c ). As 
also apparent from Fig. [2] (bottom), the Newtonian col- 
lapse produced the system with the steepest inner profile 
(7 ~ 1.7), the dMOND end-product has inner logarith- 
mic slope close to zero, while MOND collapses led to 
intermediate cases, with 7 ranging from ~ 1.2 (k = 1) 
to ~ 1.5 (k = 4). We also note that the ratio r c /rh 
(indicating the position of the knee in the density pro- 
file) increases systematically from dMOND to Newtonian 
simulations. 

The circularized projected density profiles of the end- 
products are analyzed as described in Section [3~T1 The 
best-fit Sersic indices m a , nib and m c (for projections 
along the axes a, 6, and c, respectively) are reported in 
Table 2, together with the la uncertainties corresponding 
to A\ 2 — 1; the relative uncertainties on the best-fit Ser- 
sic indices are in all cases smaller than 5 per cent and the 
average residuals between the data and the fits are typi- 
cally 0.05 < (ASB) < 0.2, where SB = -2.5 log[I(R)/I e ]. 
The fitting radial range 0.1 < R/R e < 10 is comparable 
with or larger than the typical ranges spanned by ob- 
servations (e.g., see Bertin, Ciotti & Del Principe 2002). 
In agreement with previous investigations, we found that 
the Newtonian collapse produced a system well fitted by 
the de Vaucouleurs (1948) law. MOND collapses led to 
systems with Sersic index m < 4, down to m ~ 2 in the 
case of the dMOND collapse. Figure [3] (bottom) shows 
the circularized (major-axis) projected density profiles 
for the end-products of simulations D, Ml and N together 
with their best-fit Sersic laws (m = 2.87, m = 3.20, and 
m = 4.21, respectively), and the corresponding residuals. 
Curiously, NLC06 found that low-m systems can be also 
obtained in Newtonian dissipationless collapses in the 
presence of a pre-existing dark-matter halo, with Sersic 
index value decreasing for increasing dark-to-luminous 
mass ratio. 

4.1.2. Kinematics 

We quantify the internal kinematics of the collapse 
end-products by measuring the angle-averaged radial and 
tangential components (ay and a t ) of their velocity- 
dispersion tensor, and the anisotropy parameter (3{r) = 
1 — 0.5a 2 /a 2 . These quantities are shown in Fig. [2] for 
simulations D, Ml, and N. We note that the a r pro- 
file decreases more steeply in the Newtonian than in the 
MOND end-products, while it presents a hole in the inner 
regions of the dMOND system. In addition, the dMOND 
galaxy is radially anisotropic (/3 ~ 0.4) even in the cen- 
tral regions, where models N and Ml are approximately 
isotropic (j3 ~ 0.1). All systems are strongly radially 



anisotropic for r>r^. For each model projection we 
computed the line-of-sight velocity dispersion cri os , con- 
sidering particles in a strip of width R e /A centered on 
the semi-major axis of the isophotal ellipse. The line-of- 
sight velocity-dispersion profiles (for the major-axis pro- 
jection) are plotted in the top panels of Fig. [3] The New- 
tonian profile is very steep within R c , while MOND and 
dMOND profiles are significantly flatter there. As well as 
ay, trios decreases for decreasing radius in the inner region 
of model D. The kinematical properties of M2 and M4 are 
intermediate between those of Ml and of N: overall we 
find only weak dynamical non-homology among MOND 
end-products. The empty symbols in Fig. [5] (right col- 
umn) refer to a test Newtonian simulation run with the 
FVFPS treecode (with 4 x 10 5 particles). The struc- 
tural and kinematical properties of the end-product of 
this simulation are clearly in good agreement with those 
of the end-product of simulation N (solid lines), which 
started from the same initial conditions. 

4.2. Phase-space properties of MOND collapses 

To explore the phase-space evolution of the systems 
during the collapse and the following relaxation we con- 
sider time snapshots of the particles radial velocity (v r ) 
vs. radius as in Londrillo et al. (1991). In Fig.0]we plot 
five of these diagrams for simulations D, Ml and N: each 
plot shows the phase-space coordinates of 32000 particles 
randomly extracted from the corresponding simulation, 
and, as in Fig. [TJ times are normalized to the dynamical 
time idyn (see Section [3. ip . At time t — 0.5idyn all parti- 
cles are still collapsing in simulation N, while in MOND 
simulations a minority of particles have already crossed 
the center of mass, as revealed by the vertical distribution 
of points at r ~ in panels D and Ml. At t = idyn (time 
of the peak of 2K/\W\ in the three models), sharp shells 
in phase space are present, indicating that particles are 
moving in and out collectively and phase mixing has not 
taken place yet. At t = 4td yn is already apparent that 
phase mixing is operating more efficiently in simulation 
N than in simulation Ml, while there is very little phase 
mixing in the dMOND collapse. At significantly late 
times (t = 44£dyn), when the three systems are almost 
virialized (2K/\W\ ~ 1; see Fig.[T]), phase mixing is com- 
plete in simulation N, but phase-space shells still survive 
in models Ml and D. Finally, the bottom panels show 
the phase-space diagrams at equilibrium (t = 150idyn), 
when phase mixing is completed also in the MOND and 
dMOND galaxies: note that the populated region in the 
(r,v r ) space is significantly different in MOND and in 
Newtonian gravity, consistently with the sharper decline 
of radial velocity dispersion in the Newtonian system. 

Thus, our results indicate that phase mixing is more 
effective in Newtonian gravity than in MOND 3 . It is then 

3 Ciotti, Nipoti & Londrillo (2007) found similar results in "ad 



Dissipationless collapses in MOND 



9 



-2 - 



I i i i j i i i I r i r j m r | ri r 

r b t=o.5i 



| I I I | I I I | I I I | I I I | I I 



1 1 1 j 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 




1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 

t =44t 



j i ii 1 1 ii 1 1 ii | 



1 1 1 1 1 



t = l5Ct 




.1 i i i I i i i I r i t I i i i I i 



I [ I | I M | I I I [ I I I | I [ I 

Ml i=0-5t, vn - 




[ 1 1 1 1 1 1 1 1 1 1 1 1 n r 






+H- 




-H+ 






— 




-H+ 





t=t 



<Jyn 
















1 1 [ 1 1 1 1 \i r 




H 


+H- 


-H+ 


-+++{ 



1 "T 



t=44t -: 




I I I I I I I I I I I I I I I I I H 



t = 1 50t 



— I — i — t — i — i — i — i — i — i — i — i — i — i — 



-| — i-i — ■ — | — r — i — r — | — i — i — i — | — r — i — r — | — ■ — i — ■- 




t = 4L 



| 1 1 1 | 1 1 1 | 1 1 1 | N L 


















-H+ 


-H+ 


+H- 


-H+ 


-H4 



t = 44i 



D v ■ 









+H- 


+++- 





t = 1 501 




+++ 



■J y ■ 



t i i I i i i I i i i I i I I I 



3 4.6 6 

r/r . 



3 4 , 6 8 

F /r. 



3 4,6 8 
r/r. 



Fig. 4. — Phase-space (radial-velocity vs. radius) diagrams for simulations D, Ml, and N at various times. v r is in units of v*d (simulation 
D), and u* n (simulations Ml and N; see Table 1). 



interesting to estimate in physical units the phase- mixing 
timescales of MOND systems. From Table 1 it follows 
that f* n ~ 4.7(r»/kpc) 3 / 2 (M»/10 10 M Q )- 1 / 2 Myr = 
29.8K- 3 / 4 (M»/10 10 M Q ) 1 / 4 Myr for a = 1.2 x 



10 10 ms 2 . For example, in the case of model Ml, 
adopting M* = 10 12 Af o (and r* = yjGM*/aa ~ 34kpc), 
shells in phase space are still apparent after ~ 8.3 Gyr 
(~ 44id yn )- Simulation Ml might also be interpreted 
as representing a dwarf elliptical galaxy of, say, 
M* = 10 9 M Q (and = yjGM*/ao ~ 1.1 kpc). In this 
case 44tdyn ~ 1.5 Gyr. We conclude that in some MOND 
systems substructures in phase space can survive for 
significantly long times. 

In addition to the (r,v r ) diagram, another useful di- 
agnostic to investigate phase-space properties of gravita- 

hoc" numerical simulations in which the angular force components 
were frozen to zero, so that the evolution was driven by radial forces 
only. In fact, while phase mixing is less effective both in MOND 
and in Newtonian simulations with respect to the simulations here 
reported, the phase mixing time scale in MOND is still considerably 
longer than in Newtonian gravity. 



tional systems is the energy distribution N(E) (i.e. the 
number of particles with energy per unit mass between 
E and E + dE; e.g., Binney & Tremaine 1987; Trenti & 
Bertin 2005). Independently of the force law, the energy 
per unit mass of a particle orbiting at x with speed v in a 
gravitational potential </>(x) is -E7 = i> 2 /2 + 0(x), and E is 
constant if <p is time- independent. In Newtonian gravity 
4> is usually set to zero at infinity for finite-mass systems, 
so E < for bound particles; in MOND all particles are 
bound, independently of their velocity, because (j) is con- 
fining, and all energies are admissible. This difference 
is reflected in Fig. which plots the initial (top) and 
final (bottom) differential energy distributions for simu- 
lations D, Ml, and N. Given that the particles are at rest 
at t = 0, the initial N(E) depends only on the structure 
of the gravitational potential, and is significantly differ- 
ent in the Newtonian and MOND cases. We also note 
that N(E) is basically the same in models D and Ml at 
t = 0, because model Ml is initially in dMOND regime. 
In accordance with previous studies, in the Newtonian 
case the final differential N(E) is well represented by an 
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Fig. 5. — Initial (top) and final (bottom) differential energy distribi 
and E* n /M* (models Ml and N). The energy zero points in models 
end-products have the same energy, and the highest energy particles 

exponential function over most of the populated energy 
range (Binney 1982; van Albada 1982; Ciotti 1991; Lon- 
drillo et al. 1991; NLC06). In contrast, in model D the 
final N(E) decreases for increasing energy, qualitatively 
preserving its initial shape. In the case of simulation Ml 
it is apparent a dichotomy between a Newtonian part 
at lower energies (more bound particles), where N(E) 
is exponential, and a dMOND part at higher energies, 
where the final N(E) resembles the initial one. We inter- 
pret this result as another manifestation of a less effective 
phase-space reorganization in MOND than in Newtonian 
collapses. 

4.3. Comparison with the observed scaling relations of 
elliptical galaxies 

It is not surprising that galaxy scaling relations rep- 
resent an even stronger test for MOND than for New- 
tonian gravity, due to the absence of dark matter and 
the existence of the critical acceleration ao with a uni- 
versal value in the former theory (e.g., see Milgrom 1984; 



ons. The energy per unit mass E is in units of E,j/ M, (model D), 
and Ml are such that the most bound particles of the Ml and N 
models D and Ml at t = have the same energy. 

Sanders 2000). For example, when interpreting the FP 
tilt in Newtonian gravity one can invoke a systematic 
and fine-tuned increase of the galaxy dark-to-luminous 
mass ratio with luminosity (e.g., Bender, Burstein & 
Faber 1992; Renzini &; Ciotti 1993; Ciotti, Lanzoni & 
Renzini 1996), while in MOND the tilt should be related 
to the characteristic acceleration ao- Note, however, that 
in MOND as well as in Newtonian gravity other impor- 
tant physical properties may help to explain the FP tilt, 
such as a systematic increase of radial orbital anisotropy 
with mass or a systematic structural weak homology 
(Bertin et al. 2002). Due to the relevance of the sub- 
ject, we attempt here to derive some preliminary hints. 
In particular, for the first time, we can compare with 
the scaling relations of elliptical galaxies MOND models 
produced by a formation mechanism, yet as simple as the 
dissipationless collapse. 

In this Section we consider the end-products of sim- 
ulations Ml, M2, and M4. As already discussed in 
Section [3l each of the three systems corresponds to a 
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Fig. 6. — The location of the end-products Ml (circles), M2 (squares), M4 (triangles) and M' (stars) in the planes M* — Re (a), M* - co 
(c), M, — m (d), and in the plane in which the FP is seen edge-on (b). Vertical bars account for the projection effects, while the solid 
symbols refer to average values between the two extremes. The thick solid lines represent the observed Kormendy (a), FJ (c), and FP 
(b) relations, and the thin solid lines give estimates of the corresponding scatter (Bernardi et al. 2003ab). The normalization is such that 
Rc — 4kpc for stellar mass M* = 1O 11 M0 (Shen et al. 2003). See text for the meaning of the dotted lines. 



family with constant M*/r 2 . This degeneracy is repre- 
sented by the straight dotted lines in Fig. [6^: all galax- 
ies on the same dotted line have the same k value. This 
behavior is very different from the Newtonian case, in 
which the result of a N-body simulation can be placed 
anywhere in the space R c — M*, by arbitrarily choos- 
ing M* and r„. For comparison with observations, the 
specific scaling laws represented in Fig. [5] (thick solid 
lines) are the near-infrared z*-band Kormendy relation 
R c oc M* ' 63 and FJ relation M* oc a^ 92 (Bernardi et 
al. 2003a), and the edge-on FP relation in the same band 
log i? e = A\og<jQ+B\og(M*/ Rl)+const (with ^4 = 1.49, 
B = —0.75; Bernardi et al. 2003b), under the assump- 
tion of luminosity-independent mass-to-light ratio. 

The physical properties of each model are determined 
as follows. First, for each model (identified by a value 
of k — GM*/r 2 ao) we measure the ratio i? e /?"* (see Sec- 
tion [XT]), and so we obtain R c = R c (k,M*). This func- 



tion, for fixed k and variable M», is a dotted line in 
Fig- EH- As apparent, only one pair (i? e ,M») satisfies 
the Kormendy relation for each k: in particular, we ob- 
tain that models Ml, M2, and M4 have stellar masses 
1.6 x 10 12 , 1.4 x 10 11 , and 1O 1O M , respectively (so lower 
k models correspond to higher mass systems). We are 
now in the position to obtain r* = \J GRI*/clqk and 
u* n = (ftaoGAf,) 1 / 4 , so we know the physical value of the 
projected central velocity dispersion, and we can place 
our models also in the FJ and FP planes. It is apparent 
that these two relations are not reproduced, in particu- 
lar by massive galaxies. We note that this discrepancy 
cannot be fixed even when considering the mass interval 
allowed by the scatter in the Kormendy relation (thin 
solid lines in panel a), as revealed by the dotted lines 
in Fig. [Sbc, which are just the projections of the dotted 
lines in Fig. [6^ onto the planes of the edge-on FP 4 and 

4 In Fig.JBJj the dotted lines are nearly coincident because 1) the 
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the FJ. Finally, Fig. [6ji plots the best-fit Sersic index to 
of the models as a function of M*. Observations show 
that elliptical galaxies are characterized by to values in- 
creasing with size: m > 4 for galaxies with R c > 3 kpc and 
m < 4 for those with R e < 3 kpc (e.g. Caon, Capaccioli & 
d'Onofrio 1993). Our models behave in the opposite way, 
as m decreases for increasing size (mass) of the system. 
So, while model M4 (R c ~ 1 kpc, m ~ 3.4) is consistent 
with observations, models Ml and M2 have significantly 
lower m than real ellipticals of comparable size. However, 
this finding is not a peculiarity of MOND gravity: also in 
Newtonian gravity dissipationless collapse end-products 
with to > 4 are obtained only for specific initial con- 
ditions (NLC06), while equal mass Newtonian mergings 
are able to produce high- to systems (Nipoti et al. 2003). 

So far we have compared the results of our simula- 
tions with the scaling relations of high-surface brightness 
galaxies. However, it is well known that low surface- 
brightness hot stellar systems, such as dwarf ellipticals 
and dwarf spheroidals, have larger effective radii than 
predicted by the Kormendy relation (e.g., Bender et 
al. 1992; Capaccioli, Caon & d'Onofrio 1992; Graham 
& Guzman 2003). In particular, dwarf ellipticals are 
characterized by effective surface densities comparable 
to those of the most luminous ellipticals, while their sur- 
face brightness profiles are characterized by Sersic indices 
smaller than 4 (e.g. Caon et al. 1993; Trujillo, Graham, 
& Caon 2001). Dwarf spheroidals are the lowest surface- 
density stellar systems known, and typically have expo- 
nential (to ~ 1) luminosity profiles (e.g. Mateo 1998). 
So, simulations Ml and D can be interpreted as mod- 
eling a dwarf elliptical and a dwarf spheroidal, respec- 
tively, and their end-products qualitatively reproduce the 
surface brightness profiles of the observed systems. As 
pointed out in Section 14.1.21 the velocity-dispersion pro- 
file of model D is rather flat, with a hole in the central 
regions: interestingly observations of dwarf spheroidals 
indicate that their velocity-dispersion profiles are also 
flat (e.g. Walker et al. 2006 and references therein). 

5. DISCUSSION AND CONCLUSIONS 

In this paper we studied the dissipationless collapse in 
MOND by using a new three-dimensional particle-mesh 
N-body code, which solves the MOND field equation ||3J| 
exactly. For obvious computational reasons, we did not 
attempt a complete exploration of the parameter space, 
and we just presented results of a small set of numerical 
simulations, ranging from Newtonian to dMOND sys- 
tems. The main results of the present study can be sum- 
marized as follows: 

• The intrinsic structural and kinematical prop- 
erties of the MOND collapse end-products de- 
models are almost homologous, and 2) the variable in abscissa is 
independent of k, being the FP coefficients A ~ — 2B in this case. 



pend weakly on their characteristic surface density: 
lower surface-density systems have shallower inner 
density profile, flatter velocity-dispersion profile, 
and more radially anisotropic orbital distribution 
than higher surface-density systems. 

• The projected density profiles of the MOND col- 
lapse end-products are characterized by Sersic in- 
dex to lower than 4, and decreasing for decreas- 
ing mean surface density. In particular, the end- 
product of the dMOND collapse, modeling a very 
low surface density system, is characterized by a 
Sersic index to ~ 2 and by a central hole in the 
projected velocity-dispersion profile. 

• We found impossible to satisfy simultaneously the 
observed Kormendy, Faber- Jackson and Funda- 
mental Plane relations of elliptical galaxies with the 
MOND collapse end-products, under the assump- 
tion of a luminosity independent mass-to-light ra- 
tio. In other words, this point and the two points 
above show that, in the framework of dissipation- 
less collapse, the presence of a characteristic ac- 
celeration is not sufficient to reproduce important 
observed properties of spheroids of different mass 
and surface density, such as their scaling relations 
and weak structural homology. 

• From a dynamical point of view we found that 
phase mixing is less effective (and stellar systems 
take longer to relax) in MOND than in Newtonian 
gravity. 

A natural question to ask is how the end-products 
of our simulations would be interpreted in the context 
of Newtonian gravity. Clearly, models D and N would 
represent dark-matter dominated and baryon dominated 
stellar systems, respectively. More interestingly, mod- 
els Ml, M2, and M4, once at equilibrium, would be 
characterized by a dividing radius, separating a baryon- 
dominated inner region (with accelerations higher than 
do) from a dark-matter dominated outer region (with ac- 
celerations lower than ao). This radius is ~ l.lfh for Ml, 
~ 1.8r h for M2, and ~ 2.7r h for M4. So, all these mod- 
els would show little or no dark matter in their central 
regions. Remarkably, observational data indicate that 
there is at most as much dark matter as baryonic matter 
within the effective radius of ellipticals (e.g. Bertin et 
al. 1994; Cappellari et al. 2006 and references therein). 

The conclusions above have been drawn by considering 
only simulations starting from an inhomogencous Plum- 
mer density distribution. To explore the dependence of 
these results on this specific choice, we ran also three 
simulations starting from a cold (2.K"/|W| = 0), inho- 
mogeneous and truncated density distribution p(r) — 
CM*/(rl + r 3 ), where C*" 1 = 47rln(l + r 3 /r 3 )/3, M* 
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is the total mass, and r t = 20r» is the truncation ra- 
dius. Inhomogeneities are introduced as described in 
Section [3. II Note that in the external parts the new ini- 
tial conditions are significantly flatter than a Plummcr 
sphere. The three simulations are labeled D' (dMOND), 
M' (MOND with acceleration ratio 5 k = 20) and N' 
(Newtonian). As in the case of Plummer initial condi- 
tions, also in these cases the final intrinsic and projected 
density distributions are well represented by 7-models 
and Sersic models, respectively. In analogy with model 
N, the Newtonian collapse N' produced the system with 
the steepest central density distribution (see Table 2). 
In addition, model M', when compared with the scal- 
ing laws of ellipticals (stars in Fig. [6|) , follows the same 
trend as models Ml, M2 and M4. The analysis of the 
time-evolution in phase-space of models D', M' and N' 
confirmed that mixing and relaxation processes are less 
effective in MOND than in Newtonian gravity. 

How do the presented results depend on the specific 
choice (equation [5]) of the MOND interpolating func- 
tion pi Recently, a few other interpolating functions 
have been proposed to better fit galactic rotation curves 
(Famaey & Binney 2005), and in the context of TeVeS 
(Bekenstein 2004; Zaho & Famaey 2006). It is reason- 
able to expect that the exact form of p is not critical 
in a violent dynamical process such as dissipationless 
collapse. We verified that this is actually the case, by 
running an additional MOND simulation with the same 
initial conditions and parameter k as simulation Ml, but 
adopting p(y) = y/(l + y), as proposed by Famaey & 
Binney (2005). In fact, neither in the time-evolution nor 
in the structural and kinematical properties of the end- 
products we found significant differences between the two 
simulations. This result suggests that, in the context of 
structure formation in MOND, the crucial feature is the 
presence of a characteristic acceleration separating the 
two gravity regimes, while the details of the transition 
region are unimportant. 

Though the dissipationless collapse is a very simplis- 
tic model for galaxy formation, it is expected to describe 

5 This value of re is not directly comparable with those in simu- 
lations Ml, M2 and M4, because of the different role of r* in the 



reasonably well the last phase of "monolithic-like" galaxy 
formation, in which star formation is almost completed 
during the initial phases of the collapse. The importance 
of gas dissipation in the formation of elliptical galaxies 
is very well known, going back to the seminal works of 
Rees & Ostriker (1977) and White & Rees (1978; see also 
Ciotti, Lanzoni & Volonteri 2006, and references therein, 
for a discussion of the expected impact of gas dissipation 
on the scaling laws followed by elliptical galaxies). This 
aspect has been completely neglected in our exploration, 
and we are working on an hybrid (stars plus gas) version 
of the MOND code to explore quantitatively this issue. 
We also stress that the dissipationless collapse process 
catches the essence of violent relaxation, which is cer- 
tainly relevant to the formation of spheroids even in more 
complicated scenarios, such as merging. For example, it 
is well known that in Newtonian dynamics systems with 
de Vaucouleurs profiles are produced by dissipationless 
merging of spheroids (e.g. White 1978) or disk galaxies 
(e.g. Barnes 1992) as well as by dissipationless collapses 
(van Albada 1982). Merging simulations in MOND have 
not been performed so far, and a relevant and still open 
question is how efficient merging is in MOND, in which 
the important effect of dark matter halos is missing, and 
galaxies are expected to collide at higher speed than in 
Newtonian gravity (Binney 2004; Sellwood 2004). Our 
results, indicating that relaxation takes longer in MOND 
than in Newtonian gravity, go in the direction of making 
merging time scales even longer in MOND; on the other 
hand, analytical estimates seem to indicate shorter dy- 
namical friction time-scales in MOND than in Newtonian 
gravity (Ciotti & Binney 2004). So, the next application 
of our code will be the study of galaxy merging in MOND. 
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corresponding initial distributions. 



APPENDIX 

THE VOLUME-LIMITED ENERGY-BALANCE EQUATION IN MOND 

In this Appendix we derive a useful volume- limited integral relation representing energy conservation in MOND, well 
suited to test numerical simulations. The total (ordered and random) kinetic energy per unit volume of a continuous 
distribution with density p and velocity field u is 

fc = ^(||u|| 2 + Tr4), (Al) 

where <r?- is the velocity-dispersion tensor. In the present case the energy balance equation is (e.g. Ciotti 2000) 

^ f kd 3 x = - [ pu- V0d 3 x, (A2) 
dt Jv(t) Jv(t) 
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where the integral in the r.h.s. is the work per unit time done by mechanical forces. By application of the Reynolds 
transport theorem and using the mass continuity equation we obtain 



—(k + p^ + v- [(k + P 4>)u] = P -£. 

When 4> is the MOND gravitational potential, p can be eliminated using equation |(2}, so 
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where T is defined in equation Thus, equation (|A3|) can be written as 
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By integration over a fixed control volume Vq enclosing all the system mass one obtains equation ([5]), 

THE VIRIAL TRACE W IN DEEP-MOND SYSTEMS OF FINITE MASS 



Here we prove that W = — (2/3)a/GooM| for any dMOND system of finite mass M*. Eliminating p from equation 
(fTT|) by using equation ([6]), and considering the trace of the resulting expression one finds 

W = - / 2?[^]V ■ (||V0||V0) d 3 x, (Bl) 

where we define the operator T> =< x, V >. The remarkable fact behind the proof is that the integrand above can 
be written as the divergence of a vector field, so only contributions from r — > oo are important. We will then use the 
spherically symmetric asymptotic behavior of dMOND solutions for r — > oo (BM) 



and Gauss theorem to evaluate W. 
Theorem. For a generic potential the following identity holds: 

2>[0]V ■ (||V0||V0) = V • ( X»[0]||V^||V0 - 



x||V<£|| 



Proof. From standard vector analysis (e.g. Jackson 1999) it follows that 

X>[0V • (||V0||V0) = V • (2?[0]||V0||V0) - ||V0|| < V0, VX>[0] >, 

and 

||V(/)|| < V0,VP[0] >= 



V- x||V0|| 



Identity (|B5p follows from the expansion V2?[</>] = Wcj> + V[Wcj>] as 



|V0|| 3 + HV^II < V^,2?[V^>] >= \\V<t>\\ 



V 



|V0|| 



x||V0|| 



(B2) 

(B3) 

(B4) 
(B5) 

(B6) 



3 3 

Combining equations (|B4p and (|B5|1 completes the proof of equation (|B3[) . 

We now transform the volume integral (|B1[) in a surface integral over a sphere of radius r, and we consider the limit 
for r — > oo together with the asymptotic relation (|B2[) . obtaining 



W = 



47rGao r— >oo 



4tt 3 



-gVGaoM 3 - 



(B7) 



SCALING OF THE EQUATIONS 



Given a generic density distribution p, and the mass and length units M* and r*, we define the dimensionless 
quantities x = x/r*, and p = pr\jM*. From equation ([3]) the equation of motion for a test particle can be written in 
dimensionless form as 

§ - (CD 

at 2 r% 

where 0* and t* are for the moment two unspecified scaling constants, cf> = 4>/4>* and t = t/t*, and the dimension- 
less gradient operator is V = r*V. In all of our simulations we define t» = r*/y^7, so that the scaling factor in 
equation (|C1|) is unity, while 4>* is specified case-by-case from the field equation as follows. 
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In Newtonian gravity the Poisson equation (JT|) can be written as 

'-P ■ 



V 2 <A = 4tt- 



we fix </>„ = GM*/r*, so i» = y/rl/GM* = t* n . The dMOND field equation © in dimensionless form writes 

V-(||V0||W)=47r^§^p, 



(C2) 



(C3) 



so the natural choice is = y/GM*ao, and = r*(GM*ao) 1 ' i = i*d- Finally, the MOND field equation © in 
dimensionless form is 

r i GM 

M«||V0||)V0 1 ' 



V ■ 



47T- 



where re = (p^/r^aQ. In this case, as in the Newtonian case, 0* = GM,/r„ so t* = i* n and re = GM*/r 2 ao. 



(C4) 
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